Linear scaling calculation of maximally-localized Wannier functions with atomic basis 

set 



H. J. Xiang, Zhenyu Li, W. Z. Liang, Jinlong Yang0 J. G. Hou, and Qingshi Zhu 

Hefei National Laboratory for Physical Sciences at Microscale, 
University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China 

(Dated: February 6, 2008) 

We have developed a linear scaling algorithm for calculating maximally-localized Wannier func- 
tions (MLWFs) using atomic orbital basis. An 0(N) ground state calculation is carried out to get the 
density matrix (DM). Through a projection of the DM onto atomic orbitals and a subsequent 0(N) 
orthogonalization, we obtain initial orthogonal localized orbitals. These orbitals can be maximally 
localized in linear scaling by simple Jacobi sweeps. Our 0(N) method is validated by applying it to 
water molecule and wurtzite ZnO. The linear scaling behavior of the new method is demonstrated 
by computing the MLWFs of boron nitride nanotubes. 



I. INTRODUCTION 

Wannier function^ is a powerful tool in the study of 
the chemical bonding, dielectric properties, excited elec- 
tronic states, electron transport, and many body corre- 
lations in materials. In particular, the modern theory 
of bulk polarization relates the vector sum of the cen- 
ters of the Wannier functions to the macroscopic po- 
larization of a crystalline insulator J. However, the in- 
trinsic nonuniqueness in the Wannier function definition, 
and the difficulty in defining their centers within a peri- 
odic cell calculation, limited their practical use. Fortu- 
nately, an elegant method has been recently proposed 
by Marzari and Vanderbilt to obtain a unique set of 
maximally- localized Wannier functions (MLWFs). ^ By 
transforming the occupied electronic manifold into a set 
of MLWFs, it becomes possible to obtain an enhanced 
understanding of chemical bonding properties and elec- 
tric polarization via an analysis of the MLWFs. Beside 
the above points, the MLWFs are now also being used 
as a very accurate minimal basis for a variety of algo- 
rithmic or theoretical developments, with recent appli- 
cations ranging from linear-scaling approaches^ to the 
construction of effective Hamiltonians for the study of 
baUistic transportj^ strongly-correlated electrons^^ self- 
interaction corrections, metal-insulator transitions^ and 
photonic lattices.* 

In the seminal work of Marzari and Vanderbilt, first 
a ground state calculation was carried out to obtain the 
occupied delocalized canonical orbitals, then a sequence 
of unitary transformations were performed to obtain ML- 
WFs which minimize the spread function.'^ Using the ex- 
ponential representation for the unitary transformation, 
Berghold et alS^ derived an iterative scheme to obtain 
MLWFs in large supercells of arbitrary symmetry. Also 
a simple Jacobi orbital rotation scheme was found to be 
remarkably efficient.^ A simultaneous diagonalization al- 
gorithm similar to the Jacobi diagonalization method, 
was used by Gygi et al. to compute MLWFs;'^- Zicovich- 
Wilson et al. proposed a Wannier-Boys scheme to obtain 
well localized Wannier functions in linear combination of 
atomic orbital periodic calculations. — However, all meth- 



ods mentioned above for calculating MLWFs are nearly 
O(N^) scaling (N is the number of electrons), which pro- 
hibits their applications to large systems containing hun- 
dreds or thousands of atoms. The unfavorable scaling 
comes from two steps in these methods: The conventional 
methods for getting ground state wavefunctions is 0(N'^) 
or 0(N^ InN), and the localization step in the above local- 
ization algorithms is also 0(N'^). Usually, the traditional 
ground state calculation will cost more than the local- 
ization step. However, for large systems the computing 
amount of the localization step is also time-consuming. 

In this work, we propose a simple order-N algorithm for 
effectively calculating MLWFs. The demanding ground 
state calculation is circumvented by using 0(N) den- 
sity matrix purification methods. After adopting 0(N) 
method for the ground state calculation, the conventional 
O(N^) localization step will become time-dominant for 
large systems. To obtain MLWFs in linear scaling, we 
first get initial localized orbitals from the density matrix, 
then an 0(N) localization method which uses the Jacobi 
rotation scheme is utilized to maximally localize the or- 
bitals. The linear scaling behavior of the new method is 
demonstrated by computing the MLWFs of boron nitride 
(BN) nanotubes. 

This paper is organized as follows: In Sec. we 
present our new 0(N) method for calculating MLWFs. 
In Sec. mil we describe the details of the implementa- 
tion and perform some test calculations to illustrate the 
rightness, robustness, and linear-scaling behavior of our 
methods. We discuss some possible extensions and gen- 
eralizations of our method in Sec. II VI Finally, our con- 
cluding remarks are given in Sec. 



II. THEORY 

A. Maximally-Localized Wannier Functions 

The Wannier functions are defined in terms of a uni- 
tary transformation of the occupied Bloch orbitals. How- 
ever, they are not uniquely defined, due to the arbitrary 
freedom in the phases of the Bloch orbitals. Marzari and 
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Vanderbilt^ resolve this indeterminacy by minimizing the 
total spread function of the Wannicr functions w„(r) 



(1) 



where (r)„ = (w„|r|w„), and {r^)n = {wn\r^\wn). 

Here since we aim at large systems, the F-point-only 
sampling of the Brillouin zone (BZ) is used throughout 
this work. The method of calculating MLWFs for su- 
percells of general symmetry is proposed by Silvestrelli 
et al.m^ For the sake of simplicity, considering the case 
of a simple-cubic supercell of side L, it can be proved 
that minimizing the total spread S is equivalent to the 
problem of maximizing the functional 



(2) 



where X^n = (w^m|e *"^^|w„) and similar definitions 
for Ymn and Z^n apply. The coordinate Xn of the nth 
Wannier-function center (WFC) is computed using the 
formula 



a;„ = -— Imln(w„|e ^^-''Iwn), 
zvr 



(3) 



with similar definitions for ?/„ and z„ 



B. Our 0(N) method for calculating MLWFs 

Our new 0(N) method consists of four 0(N) steps: 
first we obtain the density matrix, secondly we find 
out a set of linear independent nonorthogonal orbitals 
which span the occupied manifold, thirdly, a modi- 
fied Lowdin orthogonalization is used to orthogonalize 
these nonorthogonal orbitals, finally, the Jacobi rotation 
scheme is utilized to maximally localize the orbitals. 

In principles, any localized orbitals or density ma- 
trix based linear scaling methods can be used to ob- 
tain initial localized orbitals in linear scaling^ Here we 
use the 0(N) trace-correcting density matrix purification 
(TC2)ii method to get the density matrix since it is very 
simple, robust, and efficient. The use of some other linear 
scaling methods based on localized orbitals will be dis- 
cussed in Sec. liVI Here, the essence of the TC2 method 
is briefly outlined. In the begining, the Hamiltonian H 
represented in the orthogonal basis is normalized to an 
initial matrix pQ with all its eigenvalues mapped onto 
[0,1]: Pa = {e,naxl - H)/{e,nax - Emm), whcrc errnn and 
^max are lowest and highest eigenvalues of H, respec- 
tively. Then we correct the trace of the density matrix 
while purifying it using the following iteration: 



Pn+l{Pn 



2p„ 



Tr{p„ 
Tr{p„ 



< NJ2, 
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where Ng is the total number of electrons in a close-shell 
system. 



Given an atomic orbital, one can project out its occu- 
pied component using the density matrix operator P: 



(5) 
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where S is the overlap matrix, (j)a and cjjfi denote the 
atomic basis orbitals, P denotes the density matrix in 
the atomic orbital basis. By applying the density ma- 
trix operator on Ni, atomic basis orbitals, we can get Ni, 
localized orbitals among which only Nqcc (the num- 
ber of occupied states) localized orbitals are linear inde- 
pendent. One must select out Nocc linear independent 
localized orbitals among these localized orbitals be- 
fore performing localization steps. We have implemented 
two algorithms to achieve this goal. One of the algo- 
rithms is similar to that proposed by Maslen et 
who used Cholesky decompositions for detecting the lin- 
ear dependence. We note that in this method the total 
demanding for the construction of all the overlap matri- 
ces is 0(N) since the nonorthogonal orbitals are localized. 
Furthermore, the total computing amount for performing 
all Cholesky decompositions is almost the same as that 
for performing a sparse Cholesky decomposition with the 
matrix dimension Nocc due to the nature of the Cholesky 
decomposition. In the second algorithm, we select these 
atomic basis orbitals to be projected according to the 
physical intuition. For example, in BN nanotubes, there 
are one 2s and three 2p basis orbitals for each B and N 
atoms when using pseudopotentials. Since some electrons 
will transfer from B to N atoms, we can get Nocc linear 
independent localized orbitals by just projecting the den- 
sity matrix on all atomic basis orbitals of N atoms in de- 
spite of the large covalency in these systems. For systems 
where the bonding properties are known, the algorithm 
is found to be very efficient and the resulting orbitals are 
very sparse. In cases where the second algorithm doesn't 
apply, we will resort to the first algorithm. 

Since MLWFs are orthogonal to each other, the Nocc 
linear independent nonorthogonal localized orbitals must 
be orthogonalized. It is well known that the orthogonal- 
ized orbitals produced by the Lowdin orthogonalization 
are closest to initial nonorthogonal orbitals in the sense 
of least squares. To obtain orthogonal localized orbitals 
in linear scaling, we carry out a modified Lowdin orthog- 
onalization adopted by Stephan and Draboldii^ In this 
approach, we perform repeated first-order Lowdin itera- 
tions 



(6) 



The functions after every orthonormalization cycle have 
to be renormalized. Typically, we can obtain well orthog- 
onal localized orbitals in about five 0(N) orthonormal- 
ization cycles. 

Now we will discuss how to maximally localize these 
orthogonal orbitals to get MLWFs in linear scaling. Our 
target is maximizing ^l to obtain MLWFs in linear scal- 
ing. The key point for the successful 0(N) localization 
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step is that these orbitals are locahzed in the whole lo- 
cahzation procedure. We use the simple Jacobi rotation 
method to maximize fl since it doesn't require 0(N"^) di- 
agonalization in contrast to the unitary transformation 
method.^ This method is a traditional method in quan- 
tum chemistry for computing localized molecular orbitals 
first introduced by Edmiston and Ruedenbergr^ The ba- 
sic idea of the method is to tackle the problem of max- 
imizing Q by carrying out several Jacobi sweeps. In a 
traditional Jacobi sweep, we perform NocdNocc — l)/2 
consecutive two-by-two rotations among all pairs of or- 
bitals. The elementary step consists of a plane rotation 
where two orbitals are rotated through an angle and all 
other orbitals are fixed. In our 0(N) method, due to the 
localization of the orbitals, the computing amount of an 
orbital ratation is 0(N). Moreover, each orbital overlaps 
with only 0(1) orbitals, thus the number of the Jacobi 
rotations in a Jacobi sweep is of order N. Assuming the 
number of the Jacobi sweeps doesn't change (It is really 
the case for systems with similar characters), the total 
computing amount is 0(N). 




FIG. 1: (Color online) (a), (b), (c), and (d) show the four ML- 
WFs around a nitrogen atom of BN(5,5) nanotubes. These 
MLWFs are calculated using the supercell containing 200 
atoms with the DZP basis set. 



III. IMPLEMENTATION AND RESULTS 

A. Implementation 

Our newly developed method has been implemented 
in SIESTAfii a standard Kohn-Sham density-functional 
program using norm-conserving pseudopotentials and nu- 
merical atomic orbitals as basis sets. In SIESTA, peri- 
odic boundary conditions are employed to simulate both 



isolated and periodic systems. The details about the im- 
plementation of the TC2 method can be found in Ref. 
[18]. 



B. Validity and performance of the method 

All our calculations reported in this work are done in 
the local density approximation (LDA)>i2, Unless oth- 
erwise stated, the double-^ plus polarization functions 
(DZP) basis set is used in the calculations. First we 
calculate the MLWFs (i.e.. Boys orbitals^°) of a water 
molecule. It is well known that there are four MLWFs 
for a H2O molecule: two covalent 0-H a bonds and two 
lone-pair orbitals. The distance between the centroids 
of these four MLWFs and the position of the oxygen ion 
is 0.52, 0.52, 0.30, and 0.30 A respectively. The results 
agree well with those reported by Berghold et al..^ As a 
second check of the validity of our method, we calculate 
the piezoelectric constant of bulk ZnO. Both the Berry 
phase method^ and our new 0(N) method are used to 
calculate piezoelectric constant 633 of bulk wurtzite ZnO. 
In our 0(N) calculation, we use a 6 x 6 x 2 ZnO supercell 
since we use the F-only sampling. The results from these 
two methods agree very well: The computed values are 
1.29 and 1.30 C/m^, respectively. And both results ac- 
cord with others' result^^ (1.28 C/m^) computed through 
the density functional perturbation theory.-^ 

Then we test our method by applying it to calculate 
the MLWFs of BN(5,5) armchair nanotubes. Fig.QJshows 
four MLWFs of BN(5,5) nanotubes computed using the 
supercell containing 200 atoms. We can clearly see that 
the three MLWFs in Fig. |Il;a)-(c) are B-N a bonds. 
Among these MLWFs, the MLWF in Fig. [l^a) has ex- 
actly the same character as that shown in Fig. ^c) due 
to the mirror plane symmetry in armchair BN nanotubes. 
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FIG. 2: (Color online) Total CPU time for calculating ML- 
WFs of BN(5,5) nanotubes using the linear scaling method or 
the traditional Jacobi rotation method which doesn't take ad- 
vantage of the locahzation property of the orbitals. In case of 
the new 0(N) method, both SZ and DZP basis sets are used. 
The calculations using the traditional method are performed 
using the SZ basis set. AH calculations were carried out on a 
1.5 GHz Itanium 2 CPU workstation running RedHat Linux 
Advanced Server V2.1. 
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The other MLWF in Fig. Hid) is a tt orbital, which al- 
most centers at a nitrogen atom. It is interesting that 
in BN nanotubes MLWFs preserve a — n separation. We 
will retmn to this point later. To see the efficiency of 
our new 0(N) method, we perform a series of calcula- 
tions using the linear scaling method or the traditional 
Jacobi rotation method which doesn't take advantage of 
the localization property of the orbitals. Two different 
basis sets (single-^ (SZ), DZP) are employed in the cal- 
culations using the new 0(N) method, and only the SZ 
basis set is used for the traditional method. The CPU 
time is shown in Fig. [3 We clearly see the perfect linear 
scaling behavior of our new method. And the traditional 
method displays a nearly O(N^) scaling as expected. The 
computing saving of our method with respect to the tra- 
ditional method is dramatically large, especially when 
the size of systems exceeds 400 atoms. In addition, we 
note that the ratio between the time for calculating ML- 
WFs using DZP basis and that using SZ basis is smaller 
than the case for the ground state calculationii^ 

IV. DISCUSSION 

We have also implemented the method for obtaining 
MLWFs from the localized orbitals produced by the 0(N) 
Mauri-Ordejon (M0)2Mi or KMG^^ methods. In case 
of the MO method, the number of localized orbitals is 
equal to the number of occupied states, and all localized 
orbitals produced by the 0(N) MO method are linear 
independent. Thus the projection step is unnecessary. 
However, in the KMG energy functional, the number of 
localized orbitals is larger than the number of occupied 
states. In this case, we first get the density matrix using 
Equation (87) in the paper written by Soler et al.ii? Once 
the density matrix is available, the steps to get MLWFs 
are the same as those described in Sec. 

In our implementation, numerical localized atomic or- 
bitals are used as basis sets. Other localized basis such 
as Gaussian orbitals or real space methodsS^ can also 
be used. However, the use of non-local plane-wave basis 
which was commonly adopted in previous calculations of 
MLWFs will not result in linear scaling behavior. 

Our previous discussions mainly focus on MLWFs in 
periodic systems. We should note that our method could 
also be adopted to obtain Boy^ localized orbitals of 
isolated molecular systems with or without using peri- 



odic boundary conditions. In this case, one can directly 
minimize the spread function instead maximize the func- 
tional Vl. Besides Boys orbitals, Edmiston-Ruedenberg 
(ER)27'28 and Pipek-Mezey [VM^- localized orbitals are 
also popular among chemists. The benefit of these local- 
ized orbitals is that, unlike Boys orbitals, ER and PM 
orbitals always preserve a — tt separation. Our method 
can be used to get PM orbitals in linear scaling. How- 
ever, due to the long range character of the operator 1/r, 
our method is unable to reduce significantly the comput- 
ing amount in the calculation of ER orbitals. We notice 
that an efficient method which reduces the scaling from 
0(N5) to 0(N2)-0(N3) has been proposed by Subotnik 
et al.^ 

Although we only discuss spin-restricted systems up 
to now, combined with the spin-unrestricted linear scal- 
ing electronic structure theoryji^ our 0(N) method can 
be straightforwardly applied to insulating magnetic sys- 
tems. This method can be used to study the large multi- 
ferroic (simultaneously (ferro)magnetic and ferroelectric) 
materials. 



V. CONCLUSIONS 

To summarize, a linear scaling algorithm for calculat- 
ing MLWFs has been proposed for the first time. The nu- 
merical atomic orbital basis instead of commonly adopted 
plane wave sets is used. First we perform a linear scal- 
ing ground state calculation using the TC2 purification 
method. From the density matrix, we get the initial 
non-orthogonal localized orbitals. Through a modified 
Lowdin orthogonalization, we obtain the initial orthog- 
onal localized orbitals. Due to the localization property 
of these initial orbitals, the computing requirement of 
the subsequent Jacobi sweeps for getting MLWFs is also 
linear scaling. Our results for water molecule and bulk 
wurtzite ZnO agree well with others' results. The 0(N) 
behavior of the proposed method is clearly demonstrated 
by computing the MLWFs of BN nanotubes. Our 0(N) 
method provides a very efficient way for obtaining ML- 
WFs which have many possible applications. 
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